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The SPHC hydrodynamic code was used to sinmlate impacts of spherical aluminum 
projectiles on a single-wall aluminum plate and on a generic Whipple shield. Simulations 
were carried out in two and three dimensions. Projectile speeds ranged from 2 km/s to 10 
km/s for the single-wall runs, and from 3 km/s to 40 km/s for the Whipple shield runs. 
Spallation limit results of the single-wall simulations are compared with predictions from 
five standard penetration equations, and are shown to fall comfortably within the envelope 
of these analytical relations. Ballistic limit results of the Whipple shield simulations are 
compared with results from the AUTODYN-2D and PAM-SHOCK-3D codes presented in a 
paper at the Hypervelocity Impact Symposium 2000 and the Christiansen formulation of 
2003. 
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Nomenclature 

- Brinnell hardness of the target material 
= speed of sound 

= particle diameter marginally sufficient to produce failure 
= internal energy 
= Hugoniot value of internal energy 

- Grueneisen density factor 
= Grueneisen gamma factor 
= numerical constant 

= inter-particle smoothing length 
= mass of particle i 
= Hugoniot value of pressure 
= pressure at particle i 

- position vector for particle i 

= inter-particle separation scaled on smoothing length 
= material density 

= separation of plates of a Whipple shield 
= tensile strength of a shield plate 

* time, or, if a subscript, target 
= wall thickness in cm 

* velocity of particle i 

= relative velocity between particles i and j 
= projectile velocity 
= SPH particle mass smoothing kernel 
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I. Introduction 

T HE use of hydrodynamic simulation codes to predict behavior of engineering structures under hypervelocity 
impact has become more common in recent years. These codes are especially valuable in extending predictions 
of system behavior beyond the velocity range commonly accessible to test facilities. In order to do so in a 
believable fashion, the predictions of a given code must first be compared with those of tests, accepted analytical 
prediction methods, and other hydrodynamic codes. In particular, the effect of the formulation and code 
implementation of the equation of state (eos) must be verified before predictions at extreme impact conditions are to 
be trusted. This paper describes a set of simulations undertaken to compare results using the SPHC hydrocode 
against well-established penetration equations for single-wall impacts, and, for Whipple shields, against a semi- 
analytical expression and results from the AUTODYN-2D and PAM-SHOCK-3D codes. 

II. SPHC Code Description 

The code used far the present simulations is SPHC, a smooth particle hydrodynamics (SPH) code written in the 
C language. Each SPH particle consists of a fixed mass of material at a given location in space, together with a 
smoothing function, or kernel, that defines the particle’s extent. The half-width of the kernel is the smoothing length, 
h, and is always taken to be the same in all defined directions. In these simulations a cubic B-Spline function similar 
to a Gaussian in shape was used as the smoothing kernel, W 8 : 

Wij = K * 1.5(0.67 -? + S/2) 0<r<l 

= K * 0.25(2 - rf l<r<2 (1) 

= 0 2 <r 

Here r = | r s - rj | / h , and K = (0. 7 ntf) 1 . 

SPHC makes use of the Virtual Stress Point method of computing interparticle forces and states. In this method, 
as each particle is considered its nearest neighbors are identified and virtual stress points are created midway 
between the particle and each neighbor. Neighbor lists are normally rebuilt at each timestep using an octree 
technique. This method is faster than a simple Ar particle search, scaling as N log N. Typical initial values are about 
16 neighbors per particle in 2-D simulations and 64 in 3-D. The particles carry the mass, temperature, and velocity 
information, while the interparticle deviatoric strains, stresses, strength model parameters, and fracture info rmation 
are computed at the virtual stress points. By this method all stress forces are applied along the line of particle 
centers, rather than in directions determined by a stress tensor. Advantages of this approach are numerical stability, 
correct computation of rotations and translations, and applicability to thin shells and interfaces. It also allows a 
straightforward implementation of anisotropic strength models. 

The material density is evaluated as the sum of masses weighted by the kernel functions: 

A=Z m J W ii (2) 

j 


A “quiet start” correction begins the simulation with exactly zero pressure in all solid objects. Velocity updates are 
determined by conservation of momentum: 





Pj P\ 


\dW u 


(3) 


The pressure terms are computed from the equation of state, to which the artificial viscosity terms are added, and the 
material stress is subtracted. 

Updates to the internal energy, E, are based on conservation of energy: 


Mi 

dt 
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Here v 0 = (Vi- Vj) is the relative velocity. The P ( contains only the thermal pressure, with the artificial viscosity and 
stress terms included in the summation, divided equally between particles i and j. This is necessary to conserve 
energy exactly. The artificial viscosity used in SPHC follows the form developed by Monaghan 1 , having linear and 
quadratic terms. 

The equation of state (EOS) used in these simulations was a Mie-Grueneisen analytic form, including phase 
changes from solid to liquid and vapor, plus an ionization model. 


f , =fl,+[r 0 _+r,(— )]* p*(£-e.) <5) 

P 

The quantities P k and E* are the Hugoniot values, Yg™. is the Grueneisen coefficient, and Yi the Grueneisen density 
factor. 

The EOS computation accepts the initial material density and temperature and computes the corresponding 
pressure, internal energy, phase, and sound speed. In subsequent steps, the pressures and temperatures are computed 
from the energy as given by the energy equation. The temperature is then used to determine the phase of the 
material, and the pressure computation changes as the material melts and vaporizes. Volumetric tension in the solid 
regime is limited by the minimum pressure parameter, usually taken to be the yield strength of the material. 

At the melt temperature and above, the model approximates phase changes by a series of jumps in energy, with 
appropriate changes in the pressure law in each regime. Above the vaporization temperature die pressure is given by 
an ideal gas law. Six phase regimes are implemented: 

solid T = 0 T* u 

liquid T = T-^. -> T^ 

mol vapor T = T^ ■) Tb, (if molecular) 

atomic vapor T = T*. ■) Tio* 

ionizing T = T**, *T Mlfc 

fully ionized T>T^ 

The strength of materials model used in these simulations was a high-strain-rate-hardening model with fracture. 
This is a strain-to- fracture type of model with the fracture occurring at a specified strain and stress. The resulting 
stress/strain relation is then corrected to allow for the increased strength observed at high rates of strain. A more 
detailed discussion of the strength model is given in a companion paper in this session . Material properties needed 
in the simulations were obtained from die MatWeb.com website, http://www.matweh.com/ during the period 
October 2003 to February 2004. 


III. Single-Wall Cases 

In the single-wall cases, the projectiles were 6061-T6 aluminum, and the plates were 2024-T81 with a thickness 
of 2 mm. Projectiles impacted the plates at normal incidence. The failure criterion was the perforation or generation 
of detached spall from the back surface of the plate. 

A. Setup Details 

In all a total of thirty six individual scenarios were run. Nine were performed in three dimensions using version 
7.05 of the code on a pc-type machine from October to December of 2003. The rest of the cases were two 
dimensional, using versions 7.05 and 7.10 of SPHC on three pc-type machines from January through March of 2004. 

In making the V7.05 2-D runs, memory for up to 1,000 SPH particles was reserved, with the actual number of 
particles in a given case depending on the gridding and partitioning of the particles between the wall and projectile 
by the code. A minim um of 930 particles and a maximum of 968 particles were actually used. A “cylindrical” 
geometry was employed, so that the problem was simulated in the positive-X direction and mirrored across the X=0 
axis. The 2-mm-thick target wall was 15 or 16 SPH particles across, giving a wall particle size of 0.0133 cm or 
0.0125 on. Penetrating projectiles were from 0.045 cm to 0.106 cm, used from 14 to 60 particles each, and had 6 to 
12 particles across the diameter. This gave an SPH particle size range of 0.00750 cm to 0.01017 cm. The ratio of 
projectile particle size to wall particle size ranged from 0.60 to 0.81, so the two types of particles were of 
comparable sizes. 
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A typical run is shown in Figure 1. This is the “borderline” no-penetration case for 7 km/s. Colors are 
green=solid, light green = plastic deformation, blue = fractured, and red = melted. Note that a spall “cone” has 
begun to form in this case, but does not penetrate to the rear surface. 



For the V7.10 2-D runs memory was reserved for up to 5,000 particles. The scenario geometry was cylindrical 
as before. The actual particle counts ranged between 4,496 and 4,634. The wall was 34 or 35 particles across, giving 
a wall particle size of 0.00571 cm or 0.00588 cm. Penetrating projectiles were from 0.046 cm to 0.097 cm in 
diameter, used from 60 to 234 particles each, and had 12 to 24 particles across the diameter. This gave an SPH 
particle size range of 0.00364 cm to 0.00419 cm. The ratio of projectile particle size to wall particle size ranged 
from 0.64 to 0.73, so, again, the two types of particles were of comparable sizes. 



23456789 10 

V, km/s 

Figure 2: Single-wall spallation li mi t comparison - SPHC vs. standard penetration equations. 
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The emphasis on similar particle size is due to the dependence of simulation results on the compatibility of the 
interacting particles. If one set is much smaller than another, some of the smaller ones may “slip through the grid” 
of the larger type, producing erroneous results. Another consideration is the effect of using a cylindrical geometry, 
which is in some sense singular about the symmetry axis: particles moving toward this axis may find themselves 
constrained - by their neighbors and the impenetrability of the axis - to move along that axis. Although axially- 
traveling projectile fragments are observed in tests, the degree to which such fragments are realistically simulated in 
a cylindrically symmetric scenario is open to question. For the current set of computations, the results were 
carefully examined for axis problems, and we believe that no serious problems are present 

B. Single-Wall Results 

Minimum projectile diameters necessary to produce failure (dent) over the speed range of 2 km/s to 10 km/s were 
determined as root-mean-square values from the results of the thirty six cases. The resulting spallation limit data are 
shown in Fig. 2 and Table 1, where they are compared with predictions from five accepted penetration relations 2 : the 
Fish-Summers, Schmidt-Holsapple, Rockwell, JSC Apollo, and JSC Modified Cour-Palais penetration equations. 
The error bars on the SPHC data in the figure are + 0.006 cm, indicating the largest root-mean-square variation 
encountered among these cases. 

The SPHC results fall comfortably within the envelope given by the five standard relations, agreeing most 
closely with the JSC Modified expression fix the spallation limit: 

dZ“ = jo.01735(B//, !p r )' n (C, IV,f't w )' m (6) 


Substituting the material property values for the alloy 6061-T6 projectile and 2024-T81 plate gives: 

d J J* M = 0.19234(1^ J -0 * 316 (7) 

The SPHC results can be fitted by the power law expression: 

d s ™ c = 0.15106(V,)^ S1334 (8) 

for which the correlation coefficient is 0.99748. As seen in Fig. 1, SPHC is most conservative in the velocity range 

of 2 km/s to 3 km/s where its predictions are close to those of the Fish-Summers relation. This may be due to the 

“axial fragment” phenomenon. Above 4 km/s its predictions can be termed ‘‘moderate.” 

Table 1: Single-wall Spallation limits 

don, cm 


V, km/s 

F-S 

S-H 

Rock 

JSC A 

JSC M 

SPHC 

Fitted 

2.0 

0.113 

0.139 

0.093 

0.091 

0.124 

0.106 

0.106 

3.0 

0.081 

0.115 

0.073 

0.071 

0.096 

0.083 

0.086 

4.0 

0.064 

0.100 

0.061 

0.059 

0.080 

0.075 

0.074 

5.0 

0.053 

0.090 

0.053 

0.051 

0.070 

0.068 

0.066 

6.0 

0.046 

0.083 

0.048 

0.045 

0.062 

0.061 

0.060 

7.0 

0.040 

0.077 

0.044 

0.041 

0.056 

0.056 

0.056 

8.0 

0.036 

0.072 

0.040 

0.038 

0.052 

0.052 

0.052 

9.0 

0.033 

0.068 

0.037 

0.035 

0.048 

0.049 

0.049 

10. 

0.030 

0.065 

0.035 

0.033 

0.045 

0.045 

0.046 
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IV. Whipple Shield Cases 


In the Whipple shield cases, the shield configuration of Palmieri et al. 3 was adopted, consisting of a bumper of 
2024-T81 aluminum of thickness 1.6 mm, a standoff of 120 mm, and a backwall of 2219-T87 of thickness 3.2 mm. 
Projectiles were spheres of 2024 which impacted the bumper at normal incidence. The failure criterion was 
perforation of or detached spall from the back of the back wall. Comparisons are made with the AUTODYN-2D 
and PAM-SHOCK-3D results of Palmieri et al., and with the predictions of the generic Whipple shield penetration 
equation of Christiansen 4 . These comparisons are shown in Fig. 4 and Table 2. The SPHC results are to be 
considered as preliminary in nature, since a more detailed investigation of the materials’ behavior and the effects of 
scenario geometry is planned. 



A. Setup Details 

A total of thirteen 2-D, cylindrical symmetry scenarios were examined using version 7.07 on three pc-type 
machines during January and February of 2004. Memory for up to 6,000 SPH particles was reserved. A minimum of 
5,229 particles and a maximum of 5,531 particles were actually used. The 1.6-mm-thick bumper was 9 SPH 
particles across, giving a bumper particle size of 0.0178 cm. The 3.2-mm-thick backwall was 15 SPH particles 
across, giving a wall particle size of 0.0213 cm. The ratio of bumper particle size to backwall particle size was 0.84. 
Penetrating projectiles were from 0.30 cm to 0.66 cm in diameter, used from 132 to 572 particles each, and had 18 to 
38 particles across the diameter, which gave an SPH particle size range of 0.0160 cm to 0.0178 cm. The ratio of 
projectile particle size to bumper particle size ranged from 0.90 to 1.00; the projectile particle size to wall particle 
size ratio ranged from 0.75 to 0.84. Again, the three types of particles were of comparable sizes. 

Figure 3 shows a sequence of snapshots from the non-penetrating case velocity=7 km/s, diameter = 3 mm. few 
the Whipple Shield case. Except for the first image, which is entirely solid, the colors are green = solid, blue = 
fractured solid, and red = liquid. In the last frame, the yellow/green/red particles are a mixture of liquid and vapor 
particles At higher velocities many of the models exhibited a phase change to vapor upon collision with the lower 
shield. 

The computation can be set to put the lower layer “on hold” until late in the run, and drop the upper layers when 
they no longer interact with the downward moving material. Both of these features can be seen in the figure. In 
addition, in some runs “absorb” boundaries are used at the outer edges of the lower target to discard particles that 
have left the vicinity of the interaction regions. These techniques are used solely to save computational time and do 
not affect the final result in any way. 
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V, km/s 

Figure 4: Whipple shield ballistic limit data from SPHC, AUTODYN, PAM-SHOCK, and Christiansen. 


B. Whipple Shield Results: Medium Velocities - 3 km/s to 7 km/s 

The SPHC results are considerably more conservative at velocities between 3 km/s and 7 km/s than 
Christiansen’s expression in this range. They show a minimum in the critical projectile size of 0.30 cm at a velocity 
of 4 km/s, whereas Christiansen’s relation trends linearly upward from 3 km/s, and predicts a critical diameter at 4 
km/s of 0.40 cm. SPHC’s d^s continue to fall well below the Christiansen line for the remainder of the medium 
velocity range, being only 59 percent of Christiansen’s value at 7 km/s. 

The PAM-SHOCK data below 7 km/s only roughly bracket the critical size range. This may be due to the 
roughness of the “search grid” Palmieri et al. employed in using this code due to large run times required, lack of an 
axial symmetry option, difficulties they encountered in determining perforations reliably, and the fact that they were 
interested in examining the high velocity regime. In no case did they use PAM-SHOCK to examine projectiles 
closer together in size than 0.05 cm. PAM-SHOCK predicts no perforation by a 0.50 cm sphere at 3.3 km/s, a result 
contradicted by both the SPHC data and Christiansen’s expression. 

Neither AUTODYN nor PAM-SHOCK predicts perforation at 3 km/s for projectiles of 0.5 cm size, whereas 
SPHC’s marginal perforation diameter is 0.32 cm at this speed; Christiansen’s intermediate velocity expression 
predicts perforation by a 0.30 cm sphere. 

The crucial quantity in this velocity range is the degree of fracture of the original projectile after it encounters the 
upper shield. Experiments usually show a central fragment or fragments that will pose the main threat to the lower 
wall. The SPHC models show this behavior, but do have problems resolving the on-axis fragments due to resolution 
and numerical problems near the axis. This may lead to excess penetration. Future work with finer resolution and 3D 
models should resolve this issue. In any case, we estimate the error bars few the results in this range to be about +/- 
0.05 cm, or an uncertainty of about 1 mm in the penetration diameter. 
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C. Whipple Shield Results: High Velocities - 7 km/s to 20 km/s 

Above 7 km/s the SPHC results first approach the Christiansen curve, cross it at about 9 km/s, and then remain 
above the earlier estimates, consistently, at higher velocities. Figure 5 shows the case velocity = 14 km/s, diameter = 
5 mm at 13 jis. This is a case that should fail, according to the earlier results. Here green = solid, blue = fractured 
solid, yellow = liquid, red = vapor. The arrows indicate the direction of flow of the material at this time. The 
majority of the initially impacting material has “flashed” into a vapor or near-vapor state. This hot vapor does not 
damage the lower shield, and, in fact can interact with incoming fragments, reducing their effectiveness at 
penetrating the shield. In some tests the incoming high density material can experience a fluid instability when 
decelerated by the Iowa* density hot cloud formed from previously arriving debris, further degrading its 
effectiveness to penetrate the lower layer. The outcome is that the lower shield survives. 
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Figure 5: Case V=14 km/s, d = 5 mm, at 13 |i& 


D. Whipple Shield Results: Very High Velocities - above 20 km/s. 

For very high velocities, a simple argument holds that the details of the shield should not matter, and the -2/3 
scaling for high energy impacts should again appear. The Christiansen curve assumes this transition at 7 km/s. The 
SPHC results show some indication of approaching this limit at about 10 km/s, but the vapor mechanism described 
above could in principle shift the curve to even higher velocities. In addition, at about 25-30 km/s the upper shield 
has the capability of vaporizing the impacting projectile. SPHC has physics capable of handling high temperature 
gas, including ionization. A few tests were run with 5 mm diameter projectiles at velocities of 20-40 km/s. None 
penetrated. This means that the ballistic curve is flat out to 40 km/s. 

Since the diameter of the debris cloud when it strikes the lower plate is always equal to about the distance 
between the plates, we can easily estimate the velocity at which the momentum of a 5 mm aluminum particle will 
deliver an impulse equal to the tensile strength of the aluminum lower plate. This answer turns out to be about 50 
km/s [transit time t = S/v, implies SjS 3 = mw]. At this velocity simple momentum transfer should be able to damage 
the lower plate and the curve will begin to drop. With this scaling, d cr u depends once again on velocity to the -2/3 
power, but it varies linearly with the separation of the plates, rather than the square root dependence found at lower 
velocities. 
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Table 2: Whipple Shield Ballistic Limits 


V, km/s 



dcrit, cm 




SPHC 

Autoperf Auto no 

Pam perf Pam no Christians 

3.0 

0.32 





0.300 

3.3 



0.50 


0.50 


3.5 






0.352 

4-0 

0.30 





0.403 

4.5 






0.455 

5.0 

0.32 





0.506 

5.1 


0.80 


0.80 



5.5 



0.50 


030 

0.558 

6.0 

0.38 



0.75 


0.609 

6.5 






0.661 

7.0 

0.42 



0.75 


0.712 

7.2 



0.50 


0.50 


7.5 


0.60 

0.55 

0.60 


0.680 

8.0 

0.56 



0.60 

0.55 

0.652 

8.3 




0.75 



8.5 






0.626 

9.0 

0.60 



0.65 


0.602 

9.5 


0.56 

0.50 


0.55 

0.581 

10.0 

0.66 



0.60 


0.561 

10.5 






0.543 

11.0 

0.60 

0.50 

0.45 


0.50 

0.527 

11.5 




0.60 


0.511 

12.0 

0.60 



0.55 


0.497 

12.5 






0.484 

13.0 

0.57 

0.50 

0.45 

0.50 

0.45 

0.471 

13.5 






0.460 

14.0 

0.56 





0.449 

14.5 






0.438 

15.0 

0.56 


0.50 

0.50 


0.428 

15.5 






0.419 

16.0 






0.410 


V. Conclusion 

A preliminary set of SPHC models has been run for both single-walled and double walled meteoroid shield 
configurations. The single wall results are in good agreement with previously obtained relations. The double walled 
(Whipple) models show trends similar to earlier results, but deviate in several significant respects. First, the SPHC 
results in the intermediate velocity range (3-7 km/s) fall below the usual curve, indicating that smaller projectiles 
were penetrating than expected. Second, the SPHC results in the high velocity range (V > 7 km/s) were higher than 
the previous curve, indicating that the shields were performing better than expected and stopping larger projectiles. 
Both of these effects can be described as a general shift of the SPHC ballistic limit curve toward higher velocities. 
At least part of the high velocity effect is due to the vaporization of the lead portion of the debris cloud as it hit the 
back wall and reflected upward to interact with the remainder of the oncoming debris cloud. This is an effect not 
previously considered. Extending the effects of the vapor phase transition to higher velocities suggests that the 
ballistic limit curve may remain flat to velocities as high as 50 km/s. 
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